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1. Summary 


A general multiblock Euler solver was developed for the analysis of flow fields 
over geometrically complex configurations either in free air or in a wind tunnel. In 
this approach, the external space around a complex configuration was divided into 
* a number of topologically simple blocks, so that surface-fitted grids and an 
efficient flow solution algorithm could be easily applied in each block. The 
computational grid in each block is generated using a combination of algebraic and 
elliptic methods. A grid-generation/flow-solver interface program was developed 
to facilitate the establishment of block-to-block relations and the boundary 
conditions for each block. The flow solver utilizes a finite-volume formulation and 
an explicit time-stepping scheme to solve the Euler equations. A multiblock 
version of the multigrid method was developed to accelerate the convergence of 
the calculations. The generality of the method was demonstrated through the 
analysis of two complex configurations at various flow conditions. Results were 
compared to available test data. Two accompanying volumes, user manuals for the 
preparation of multi-block grids (V ol. II) and for the Euler flow solver (V ol. Ill), 
provide information on input data format and program execution. 
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2. Introduction 


The increases in fuel cost over the last two decades has spurred the 
development of fuel-efficient propulsive devices for transport aircraft. Studies 
(Ref. 1 , 2) have shown that high speed propfans have the potential of reducing the 
fuel consumption of conventional subsonic transports. More recently, both fuel 
efficient advanced turbofan and very-high-bypass superfan engines have been 
developed for a new generation of advanced transport aircraft. 

These propulsion systems offer significant reduction in specific fuel 
consumption. It is essential, however, to integrate the propulsion system properly 
with the airframe if the fuel-efficiency advantages are to be realized and optimized 
in actual flight. NASA Langley Research Center has an on-going research 
program on the installation of advanced propulsion systems on transport aircraft 
with supercritical wings (Ref. 3). Complementary to that effort, the present work 
is concentrated on the development of a computational fluid dynamics (CFD) tool 
for studying the effects of advanced propulsion systems on airplane lift and drag 
characteristics. 

In developing a CFD tool for engine-airframe integration, one must address the 
following issues: complex geometry; nonlinear effects associated with transonic 
flows, rotational flows and flows with non-uniform total pressure and temperature 
(caused by propfan or turbofan engines); and viscous effects on the configuration 
surfaces. Ideally, a complex flow such as this can only be resolved with the 
Navier-Stokes equations. Unfortunately, the Navier-Stokes formulation requires 
very fine grids to resolve the viscous layers, which, in turn, requires large amounts 
of computing resources in order to obtain useful results for the analysis of a 
complex airplane configuration. An alternative is to neglect all viscous terms and 
reduce the governing equations to the inviscid, Euler equations, so that they can be 
solved effectively with the available computing resources. 

In this report, an Euler-based computational method for the study of engine- 
airframe integration problems will be presented. Two major tasks were involved 
in the development of this method; the generation of computational grids; and the 
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development of an appropriate flow solver. 

In the area of grid generation, extensive research can be found in the literature. 
These include the unstructured grid approach (Ref. 4, 5), and the structured grid 
approach (Ref. 6, 7). The unstructured grid approach is an ideal choice for 
complex configuration analysis. However, the requirement for large central 
memory and the relatively low execution efficiency, presently limits the practical 
applications of the unstructured grid method (Ref. 8-11). The structured grid 
approach together with a multiblock method to handle complex geometries (Ref. 
12-17), offers both geometry flexibility and execution efficiency for large 
problems. This approach was chosen as the framework for the present study. 

The basic philosophy behind the multiblock grid generation process is 
discussed in Section 3. A grid generation toolkit, program BCON, was developed 
to facilitate this process. BCON is documented in details in volume II (Ref. 18) of 
this report. 

The Euler-based flow solver was developed for general configuration analysis. 
Because of the multi-block concept, this General Multi-Block Euler (GMBE) 
solver is both configuration independent and grid topology independent. The 
complete flow field, exterior to the configuration, is divided into a number of 
topologically simple regions or blocks. Each block consists of six faces in the 
computational plane. Along each face of the block, boundary conditions, together 
with block-to-block relationship; are specified. Every block has similar data 
structure, so that the same flow solver algorithm can be used throughout the entire 
flow field. 

Two new features are introduced in the present multiblock flow solver: 1) 
generalized block- to-block connectivity which allows each block an arbitrary 
number of neighbors, and 2) general boundary conditions, to allow each face an 
arbitrary combination of boundary conditions. These new features relieve the 
blocking restrictions encountered in most existing multiblock methods, where each 
face of the block can only have one neighbor and one type of boundary condition. 
These features result in a significant reduction in blocks and data communication 
requirements. As a result, there is lower block boundary overhead and more 
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effective vectorization in the flow solver. These two features can also be very 
effective for simpler configurations using very few (2-10) blocks (Ref. 19). 

Computed results for two test cases are presented to demonstrate the generality 
of the method. The first case is a wing-mounted propfan configuration in a wind 
tunnel. The second case is a generic twin-engined, NASA low-wing transport 
configuration (Ref. 3). Results from both are compared with test data. Volume HI 
(Ref. 20) is the user manual for the Euler solver. 


3. Multiblock Grid Generation 

The generation of a surface-fitted grid over a geometrically complex 
configuration, in a manner that is non-configuration specific, is an essential 
element in the present general flow solver development. In general, a single grid 
cannot provide adequate topological resolution to cover a complex shape. 
Consequently, the flow field needs to be divided into a number of simple regions, 
as illustrated in Figure 1, such that within each region, smooth and well distributed 
grids with a desirable grid topology can be generated. Each block contains six 
faces in the computational space. Some faces may have no neighbor, such as the 
faces on the far-field boundary or on the configuration surface. Others may have 
multiple neighbors, such as those on the interface boundaries. The configuration 
surface must lie on block-faces. This is necessary for a non-configuration specific 
multi-block algorithm. Along each block-face, a combination of boundary 
conditions may be used. Presently, there is a restriction that the common faces of 
adjoining blocks have identical grid points. 

The multiblock grid generation process is divided into four steps: 

1. Blocking of the Complete flow field - . > 

The first step in the grid generation procedure is to divide the complete flow 
field into a number of topologically simpler blocks. Complex geometry, such as a 
complete airplane, may require many blocks (currently, of the order of 30). 
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2. Surface Grid Generation on Block Faces 


In this step, the grids on each block face must be generated. If the block face is 
part of a configuration surface, the grid must be interpolated in the parametric 
space of the original geometry definition to ensure that the grid lies on the 
configuration surface. All other block faces will be either an interblock boundaries 
or the far-field boundaries or inlet or exhaust boundaries or propeller disk 
boundaries. The surface grid for these faces may be generated either by an 
algebraic method or by an elliptic grid generation method. 

3. Volume-Grid Generation 

After the surface grids of each block have been properly prepared, the EAGLE 
(Ref. 9) grid generation package can be used for field-grid generation. This 
package uses both algebraic/transfinite interpolation and elliptic grid-generation 
methods for volume-grid generation. 

4. Block-to-Block Relation and Boundary Condition Definition 

The block-to-block relation file is required by the flow solver. This file 
establishes the flow field communications between adjacent blocks. In addition, a 
block-boundary condition file in needed by the flow solver in order to apply the 
proper boundary conditions along each block face. The block-to-block relations 
are established such that every common block boundary is assigned a unique 
record number, together with the two dimensional grid index associated with that 
boundary. The record number of its direct neighbors are also identified. This 
information enables the Euler solver to access appropriate data from neighboring 
blocks for the convection and dissipation flux calculations. The block-boundary 
condition is assigned in patches. Each patch represents a unique boundary 
condition with a given 2D index. Each block face can have any combination of the 
following types of boundary conditions: 1) interface, 2) solid surface, 3) inlet, 4) 
exhaust, 5) far-field 6) centerline, 7) mirror image, and 8) actuator disk. 

BCON is used to prepare the geometry data and job deck inputs to run the 
EAGLE code (Ref. 6) for volume grid generation. BCON will also generate the 
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block-to-block relation and block-boundary condition files for the Euler solver. In 
addition, BCON can be useful for data checking to ensure that all information is 
correct, before it is passed to the flow solver. 

4. General Euler Solver Algorithm - 

The Euler solver was derived from Jameson’s cell-centered finite-volume 
approach (Ref. 21) where the unsteady Euler equations are solved using a multi- 
stage Runge-Kutta time-integration method; see Appendix A. The main objective 
of this work was to develop a flow solver that can accept any combination of grid 
topologies for the analysis of complex flows over arbitrary configurations. 
Jameson’s solver uses a C-grid-topology which is adequate for the analysis of 
wing/body flow. To achieve the objective of treating arbitrary geometries, the flow 
solver must be non-configuration specific and non-grid-topology specific. This 
objective required the adaptation of Jameson’s flow solver to a multiblock grid. In 
the multiblock version of the Euler flow solver, boundary condition surfaces are 
restricted to the block faces (Sec. 3). Then the interior of every block would 
contain only regular flow field grid points, such that the flow field of each block 
can be solved using the same flow solver. Requirements for the adaptation of 
Jameson’s flow solver to multiblock grid include: development of a multiblock- 
compatible version of the multigrid method; algorithm improvement for treating 
irregular and singular grid points in the flow field; and implementation of boundary 
conditions in the multiblock environment. A flow solver that addressing all these 
issues is necessarily demanding on computing resources. Thus efficient 
management of computer memory is an important aspect of GMBE. 

The multigrid strategy employed in the present study is a simple V-cycle 
multigrid applied within each individual block of the flow field. Two main reasons 
for using V-cycle multigrid in a block-by-block manner are: 1) to minimize data 
transfer from main memory to the secondary memory (Ref. 22) and 2) to 
significantly reduce the programming complexity for the present multiblock solver. 
However, this restricts the multigrid calculation to one block at a time. The 
residual collection and the correction interpolation at the grid cell near the block 
face can not be done in the same way as the single block approach. This can affect 
the overall convergence. It may also cause numerical oscillations for the analysis 
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of supercritical wings with rooftop-type pressure distributions (Ref. 12). A simple 
way to circumvent this problem is to turn-off the multigrid calculation after the 
solution converges to a certain level of accuracy. Iterative computations will then 
be continued in the single grid mode until desired convergence is achieved. In 
dealing with complex geometries such as the wing/body/under-wing nacelle/pylon 
configuration shown in Fig. 2, it can be very time-consuming to ensure that every 
block of a multiblock grid does not contain an odd number of cells in some given 
direction within the computational space. An adaptive multigrid strategy has been 
developed which recognizes blocks which have an odd number of grid cells in any 
given direction. A block containing an odd number of grid cells, in any direction, 
will be solved by single-grid calculations. Different multigrid levels in other 
blocks are employed based on the local grid cell distributions. In order for the 
adaptive multigrid strategy to be effective, the non-multigrid blocks should be 
small in number and size; most of the flow field should be solved using multigrid. 

Jameson’s cell-centered finite-volume scheme (Ref. 21) in the present Euler 
solver is accurate and reliable provided that the grids in the flow field are smooth 
and well distributed. However, irregular and singular grids in the flow field can 
affect the accuracy and the reliability. Examples are the collapsed edge in front of 
the wing leading edge near the wing tip (Fig. 3a), and the fictitious comer along 
the fan cowl surface (Fig. 3b). Special techniques must be applied in these regions 
in order to obtain an accurate and reliable solution. There are two main concerns 
here: 1) to ensure flux balance; and 2) to minimize the sensitivity of the dissipation 
terms to grid irregularities such as high aspect ratio cells and fictitious comers. 
The convective flux and the dissipative flux along these grid lines must be fully 
conservative. Flux terms along each face of the cell are computed and used 
consistently for the flux balance of the two cells adjacent to the face. Treatment of 
zero volume cells near the collapsed edge is given in Appendix A. Another 
improvement that has been added to the flow solver is the new dissipation 
formulation based on a spectral radius scaling (refs. 23-25). The new formulation 
improves the consistency of the solution by decreasing sensitivity to the dissipation 
parameters. The improvement over conventional dissipation formulations (Ref. 
26) is especially pronounced on grids with a large aspect ratio or with nonsmooth 
mesh distributions (Ref. 27). Further details can be found in Appendix B. 
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On each face, the calculation process is dependent upon the type of boundary 
condition specified. If it is an interblock boundary surface, conservation of mass, 
momentum, energy, and dissipation flux are ensured across the block boundary. 
The fourth-order dissipation requires two layers of data from every adjacent block 
in order to formulate a flux term that is identical to the one used in single-block 
calculations. 

Characteristics boundary conditions (Ref. 28) have been implemented along the 
far-field boundaries. This ensures proper propagation of information in and out of 
the boundaries of the flow domain. This is crucial for numerical accuracy and 
stability, especially for cases, involving powered nacelles or propellers, with non- 
uniform total temperature. The tangency boundary condition is implemented on 
the configuration surface, where the normal momentum equation is used to 
compute the surface pressure. 

Propeller power loading is simulated by an actuator disk where the 
distributions of total pressure, total temperature, and swirl are prescribed. 
Another option for simulating propeller power effects is to prescribe the thrust, 
normal force, and side force on the propeller disk and the work added to the flow 
by the propeller. The flow variables downstream of the disk are related to their 
upstream values through the continuity, momentum, and energy equations. The 
major advantage of this method is that the effects of angle of attack, as well as the 
influence of side flow, can easily be simulated through the input of normal force 
and side force distributions on the propeller disk. This method is described in 
more detail in References 29 and 30. 

Boundary condition methods for turbofan engine inlet and exhaust are given in 
Reference 31. The nacelle center line, as well as any grid line in the flow field 
representing a collapsed grid plane, is treated as solid boundary with zero face 
area. Also sufficient ghost cells are employed in the mirror-image boundary 
conditions for producing exact symmetry solutions. Further details on boundary 
condition implementations are given in Appendix C. 

The efficient management of directly addressable central memory is a necessity 
in large scale CFD simulation because central memory is a scarce resource in a 
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multi-processor supercomputer. Therefore, for a general flow field simulation, it 
would be very useful to adjust the memory requirement to each grid for a specific 
CFD application. To facilitate this process, a pre-processor to the general 
multiblock Euler code was written to manage the central memory requirement for 
the Euler calculations. Grid information is read in and processed by this pre- 
processor code to ensure optimum block-storage area in the memory for storing the 
flow field and geometry information of each block. This information includes the 
flow solution vector, grid coordinator vector and cell-volume vector. The 
information is then moved to a working memory location, in a block-by-block 
manner, to allow the flow solver to update the solution. The updated solution is 
then stored back in the block-storage area. The working memory is sufficiently 
large to hold the data for any one block including dissipation vectors, and other 
temporary arrays required for the Euler code multigrid calculations. 


5. Results and Discussions of Test Cases 

5.1 Case 1: Wing-Mounted Propfan in Wind Tunnel 

A wing-mounted propfan configuration in a wind tunnel. Fig. 4, was analyzed at 
Moo = 0.167, a = 0.0 deg. Solid wind tunnel wall boundary conditions are 
included in the calculation to account for blockage effects. The propeller power 
input (in terms of total pressure, total temperature, and swirl) is prescribed on the 
propeller disk for a low power setting. As the general multiblock Euler code is 
neither configuration nor grid-topology specific, it allows for the selection of a grid 
system that is optimum for a specific application. For this test case, a cylindrical 
type of grid is used to provide good grid resolution on the propeller disk and 
nacelle surface. Fig. 5 shows the configuration and several slices of the grid used 
in the analysis. The grid slices in Fig. 5c and 5d show the dense grid near the solid 
body boundary on the nacelle and wing, as well as near the tip of the propeller 
disk. The complete flow field is divided into two zones by a cylindrical-type grid 
surface. At the propeller disk station, the cross-section of this surface coincides 
with the tip of the propeller disk. The inner grid zone between this surface and the 
configuration surface is divided into nine blocks. The outer grid zone is also 
divided into nine blocks. The field grid consists of 860,000 points. 
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The isobars on the configuration surface for the power-on case are shown in 
Fig. 6. The effects of propeller power are depicted in Fig. 7 by comparing the 
propeller power-on (solid line) with power-off (dashed line) results. Wing section 
lift has been decreased on the downwash side due to the effects of swirl, Fig. 7(a). 
This is accompanied by an increase of section lift on the upwash side. Fig. 7(b). 
The computed section pressure coefficient (Cp) distributions showed excellent 
agreement with test data for both the downwash. Fig. 8(a), and upwash Fig. 8(b), 
locations. 


5.2 Case 2: Low Wing Turbofan Transport 

To further demonstrate the code’s capability to handle complex geometries, a 
NASA low wing transport configuration consisting of wing/body/under-wing 
nacelle/pylon was analyzed. This model was tested in the NASA Langley 16-foot 
Transonic Tunnel (Ref. 3) for the investigation of nacelle/pylon installation on 
wing/body. Detailed information about the model can be found in Ref. 3. The 
surface grid is given in Fig. 2. The entire flow field consisted of 29 blocks with a 
total of 1,140,000 grid points. The geometric complexity required the use of a 
composite grid system for proper resolution of the flow field around the 
engine/nacelle and also near the wing. The composite grid system included an FI- 
type grid for the external flow field and a cylindrical grid for the nacelle inlet and 
exhaust flows (Fig. 9). To allow for the study of engine power effects, the detailed 
gridding around the nacelle was constructed with a grid plane on the inlet fan-face, 
a grid plane on the exhaust face of the primary (or core) exhaust flow and another 
grid plane on the exhaust face of the fan exhaust flow. To simulate a flow-through 
condition, the total pressure and total temperature at both exhaust planes were set 
at freestream values. The velocity at the fan inlet face was arbitrarily set at 0.6 
times the freestream velocity to account for inlet blockage. Fig. 10 shows the 
- isobars on the configuration surface for the flow-through case when Me = 0.77 and 
a = 0.5 degree. It can be seen that detailed flow resolutions in the nacelle -pylon 
region are obtained with the present multiblock Euler solver. Figure 1 1 compares 
the computed wing surface pressures at two span stations inboard of the nacelle 
with test data. The differences could have arisen because the viscous effects and 
the aeroelastic effects were not accounted for in the present analysis. 
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This configuration was also analyzed with engine power-on to demonstrate the 
code’s capability to simulate power-on conditions. The total pressure Po 
(normalized by the freestream static pressure) for the core exit was set at 1 .92 and 
at 1.60 for the fan exit. The inlet velocity at the fan face was set at 0.6 times the 
■ freestream velocity. Figure 12 compares the Cp distribution, for the engine 
power-on and power-off cases, along the outboard intersection line between the 
pylon and the primary cowl. The difference in the pylon surface pressure 
distributions are discemable even at this low power setting. ....... 

Finally, the same configuration was analyzed at an off-design condition (M oc ,= 
0.8, and a = 0.474 deg). The results are compared with test data in Fig. 13. The 
discrepancy between experimental and computed shock locations are more 
pronounced because of the stronger viscous effects at this higher transonic Mach 
number. 


6. Concluding Remarks 

A grid generation and block-to-block boundary condition pre-processor 
(BCON) and a general multiblock Euler (GMBE) solver were developed. The 
generality of the program is demonstrated through the analysis of two complex 
configurations at different flow conditions. The analysis results agree well with 
test data provided that the viscous effects and aeroelastic effects are not 
predominant. When used within these conditions, the code is a useful and 
effective analysis tool for engine/airframe integration studies. Accompanying 
volumes are user manuals for BCON (Vol. II) and GMBE (Vol. HI) 
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Appendix A 

Numerical Method for the Euler Solver 


The basic numerical scheme follows Jameson’s finite volume formulation (Ref. 
26). In this formulation, the computational domain is divided into a large number 
of computation cells each treated as a control volume for the flux balance of the 
governing conservation laws. For the Euler equations 


3w/3t + 3f/3x + 3g/3y + 3h/3z = 0 (Al) 
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Here, the conventional notations p, u, v, w, p, E are used for density, three 
components of velocity, pressure, and the total energy. The total enthalpy H is 
defined as 


H = E + p/p = yp/(y -l)p + (u 2 + v 2 + w 2 )/2 


(A2) 


with y representing the ratio of specific heat. Equations (Al, A2) are augmented 
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by the perfect gas law as the thermal equation of state, and the caloric equation of 
state for the perfect gas. In the governing equations, the density and the pressure 
are normalized with respect to their freestream values. The total energy and the 
total enthalpy are normalized with respect to the ratio between the freestream 
pressure and freestream density poo/p<», and the velocity is normalized with respect 
to the square root of poo/p«>. 

The unsteady Euler equations are solved by the multistage Runge-Kutta time- 
stepping scheme which for an n stage may be written (Ref. 23) for the time step n 
to n+1 as: 


w (0) = w n 

w 0) _ w (°) - ai AtR(0) 

(A3) 

w (m-l) _ w (0) _ (x m _ 1 AtR^" 1-2 -* 

w n+1 = w*-™) = w^ - oc m AtR^ m-1 ' ) 
with R (m) = C(w) (m) + D(w) (m) 


Here, C(w)^ represents convective flux terms, and D(w/ m ^ represents artificial 
dissipation terms. The convective terms {f } , {g} and (h) from Equation (Al) are 
discretized as given in Reference 26. The dissipation terms are added to control 
numerical oscillations and to capture the shocks. Details on the dissipation 
formulation can be found in Appendix B. Central difference approximations are 
used for both convective and dissipative term calculations. 

The implicit residual smoothing (IRS) method of Ref. 32 was employed to 
extend the stability range of the Runge-Kutta explicit time marching scheme. IRS 
also significantly enhances the smoothing characteristics of the time marching 
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scheme for use with multigrid. This version of IRS provides superior 
performance, particularly for very fine grid calculations. 

In addition, for flow with uniform enthalpy, enthalpy damping procedure (Ref. 
26) is used to accelerate the convergence of Euler solutions. It modifies the Euler 
equations by adding forcing terms proportional to the difference between the local 
and free stream values of total enthalpy. The formulation of spatial discretization 
ensures that constant total enthalpy condition is consistent with the steady-state 
solution of the difference equations. In so doing, the steady-state solution is not 
altered by these additional terms. 

Many grid surfaces converged into a collapsed edge as shown in Fig. 3a. 
Extending of such coinciding surfaces beyond the collapsed edge would result in 
zero volume cells. These cells must be eliminated from the block grid definition 
because flow physics requires that the adjacent blocks on both sides of the zero 
thickness surface access information from each other directly. This was acheived 
by specifying a special relationship between the adjacent blocks. The information 
is given in the block-to-block relationship file that was created by BCON (V ol. II) 
during the process of multiblock grid generation. 
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Appendix B 

Artificial Dissipation 


The artificial dissipation model used here is based on the model devised by 
Jameson and Baker (Ref. 21) for 3-D Euler equations. The original model was 
modified for scaling in the coordinate directions as suggested by Martinalli (Ref. 
25) and Swanson and Turkel (Ref. 24). These concepts are outlined below. 

A spectral radius scaling on the dissipation terms is used to minimize the 
amount of artificial dissipations. The dissipation term D(w) is expressed as 


D(w) = D(w)j + D(w)j + D(w)k 


(Bl) 


where the indices i, j, and k denote contributions from the rj, and £ directions 
respectively. The formulas for the contribution from ^-direction are discussed 
here. The corresponding formulations in the two other directions can be obtained 
through appropriate permutation of the indices. 


D(w)j = 

di+l/2,j,k ■ di-i/2,j,k 


(B2) 


These dissipation terms are a blend of second and fourth order differences with 
coefficients adapted to the flow. 


<W.j.k = Bi i+ i /2 ,j.k(E <2> - e (4 >8 5 2 >(W i+! .j,k - Wi,j.k) 


(B3) 


where the scaling factor Bii + ]/2,j,k = ^i+i/2,j,k is the spectral radius in the £- 
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direction at face (i+1/2, j, k). Bi is first computed at cell center (i, j, k) 


Bii,j,k - M.j.k - (Any 


i,j.k 


+ a i,j,k)(Si-l/2,j,k + S i+l/2,j,k)/2 


(B4) 


where ajjk is the speed of sound at the cell center and Sj+i/2 t j,k is the cell face 
area. Bi at the cell face (i+1/2, j, k) is then computed by averaging the two values 
at the adjacent cell centers, i. e.. 


Bii+i/2,j,k “ (Bij.j.k + Bij+ij t k)/2 (B5) 

In Equation (B3) e^ 2) and e^ 4 -* are the adaptive coefficients and the 5^ 2 is the 
second difference operator in the ^-direction. The coefficient e (4) controls fourth 
order dissipation which provides the background dissipation in the smooth part of 
the flow to suppress numerical oscillations. Near shocks it is reduced to zero. The 
coefficient e (2) controls second order dissipation which is used to capture the 
shocks. The term e (2) is made proportional to the normalized second difference in 
the pressure 


e (2) = k ( 2) max(vijj f k, vi j+1>j>k ) , 

^M.j.k = KPi+l,j,k " ^Pjj^k + Pi-1, j,k)/(Pi+l,j,k + 2pij t k + Pi-1, j,k)l > 
e (4) = max[0, (K^ - e®)] , 


(B6) 

(B7) 

(B8) 


where typical values of the constant K (2) and K (4) are 1/4 and 1/256, respectively. 
This isotropic dissipation formulation leads to good shock capturing and is 
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adequate for computational grids of moderately high aspect ratio (of the order 10). 
However, for grids containing very high aspect ratio grid cells (e. g., 500), the 
dissipation coefficients may become imbalanced in different coordinate directions. 
In these cases, the artificial dissipation should be anisotropic in the three 
coordinate directions. According to the modification developed by Martinelli (Ref. 
25), the scaling factors of Equation (B4) are adjusted as follows: 


Bii,j,k = <t>(ri)Wi,j,k ® 9 ) 

where 

ri = (^ji,j,k ■*" ^q,j,k)/^i,j,k 
and 

O(ri) = 1 + ri a where 0 < a < 1 


a = 2/3 is recommended by Martinelli. 

The pressure sensor term based on the normalized second difference of pressure 
in Equation (B7) has some interesting characteristics. Across a shock wave, this 
term will be large since a large pressure variation can lead to large second pressure 
difference. In the smooth part of the flow, this term will be insignificant. For 
some low speed flow calculations where the shock wave does not occur, this 
pressure sensing term may become significant near the leading edge, where the 
flow is undergoing rapid expansion. For this type of flow, the pressure sensor of 
Equation (B7) is replaced by a sensor using an entropy function f defined as 


f = (P/P) Y 
where 
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Vij.j.k — l(fi+l,j,k " ^i,j,k + fi-l,j,k)/(fi+l,j,k + ^fjjk + fi-l,j,k)l 


(B9) 


The value of this sensor term will be small in regions of flow where the entropy 
function f (or entropy) is uniform or a slow varying function. 

In the coarse grid a simplified dissipation model is used by replacing Equation 
(B3) with 

d i+ i/2,j,k = K' 0) Bi i+1/2J , k (W 1+1 ,j, k - W,j. k ) (BIO) 

where typical value of the constant is 1/128. 
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Appendix C 

Boundary Condition Implementations 


The boundary conditions discussed in this Appendix include: 1) Solid-surface 
boundary condition, 2) Far-field boundary condition, 3) Nacelle-inlet boundary 
condition, and 4) nacelle exhaust plane/propeller disk boundary conditions. 

1). Solid-surface boundary condition 

Along the configuration surface, the tangency condition is enforced by setting 
q * n = 0 where q is the velocity vector and n is the unit outward normal vector. 
For the cell-centered scheme used in the present work, the pressure on the 
configuration is also needed to evaluate the convective flux terms in the 
momentum and energy equations. Simple extrapolation from adjacent cell centers 
(see Vol. in pages 8 and 16) can be used. Optionally, to improve numerical 
accuracy, a normal momentum equation that relates the normal derivative of 
pressure to the streamwise and spanwise derivatives can be used along a solid 
surface in a plane of constant j the relation becomes 


s^ ap/a^ - pUiUjacsjj)/^ = o 


(Cl) 


with Sy = h a^i/axj , the projected area of cell face-i in the xj direction, where xj 
= x , x 2 = y , x 3 = z , and = £ , £ 2 = J] , £ 3 = C , h = laxj/aSjl , the cell 
volume, and Uj = SyUj , the contravariant velocity component in the ^ -direction. 

2). Far-field boundary condition 

Because the calculation is performed in a finite region, correct implementation 
of the far-field boundary condition is essential for numerical accuracy and stability. 
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The characteristic boundary condition based on locally one-dimensional Euler 
equations is used. Along each face of the far-held boundary cell, the velocity 
component normal to the far-held boundary is computed and checked. 

2(a) For a subsonic freestream condition, if it is an inflow boundary, there are 
four incoming characteristics requiring the specification of four quantities: the 
tangential velocity; the enthalpy; the entropy; and one Riemann variable evaluated 
from upstream values. This information together with another Riemann variable 
extracted from the interior flow held are used to determine all flow variables, i.e., 


qt — q^x, > s — Soo , H Hoo , 

q„ - 2a/Cy- 1) = qne - 2ae/(y- D , (C2) 

qn + 2a/(Y- 1) = qn_ + 2a«/(y- 1) . 


where s denotes the entropy, a denotes the sonic speed, q the velocity, subscript e 
representing values extrapolated from the interior cells adjacent to the boundary, 
subscript °° denoting freestream values, subscript t implying tangential, and 
subscript n implying normal components. 

2(b) For the subsonic outflow boundary, specification of one quantity is 
required for the incoming characteristic. One Riemann variable based on 
freestream values together with four extrapolated quantities are used to calculate 
the flow variables, i.e.. 


qt — q^ » s — s e , H , 

q„ - 2a/(y- 1) = q^ - 2a co /(Y- 1) , (C3) 

q„ + 2a/(y- 1) = + 2a e /(y - 1) . 
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This approximation assumes a uniform total temperature distributions down 
stream in the far-field. This assumption is incorrect for powered nacelle. 
Therefore, for powered nacelles an alternative downstream far-field boundary 
condition is specified. The conservative variables p, pu, pv, pw, and pE are 
extrapolated from their upstream values, and the pressure is set to the freestream 
value. 

2(c) For a supersonic inflow boundary, all characteristics are incoming. Hence 
all variables are assigned to their freestream values. 

2(d) For a supersonic outflow boundary, there are no incoming characteristics. 
Hence, all flow variables are extrapolated from the interior flow solutions. 

3). Nacelle-inlet boundary conditions 

These boundary conditions are described in more detail in Ref. 31. The nacelle 
inlet face is treated as an outflow boundary. This is possible because the interior of 
the nacelle is not modeled. Three choices of boundary conditions exist: 

a). Mass flux boundary condition 

If the capture stream tube area of the nacelle, i.e., A„> is specified, the mass 
flux at the inlet face is given by 


(PQn)f (PooAoo(Joo)/Af 


(C4) 


with Af representing inlet face area. This is equivalent to specifying an inflow 
area ratio. 

b). Velocity boundary condition 

For a given captured stream tube area Aoo , the mass flux boundary condition 
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can be converted into a velocity boundary condition using the isentropic relation to 
relate density to the local velocity, i.e., 


qn/q« = (A_/A f )[l - (y-l)(M„) 2 {(q n /q„) 2 - \\Hr V<1 ~ V> (C5) 


A simple Newton’s iteration procedure is written to solve the above nonlinear 
equation for q^ When shock waves are presented in the inlet, the accuracy of this 
equation deteriorates as the shock waves get stronger. A well designed inlet for 
subsonic transport, however, should not have strong shock waves in the inlet. Test 
Case 2 in Section 5.2, using qn = 0.6qoo, is equivalent to using an inflow area ratio 

Aoo/Af = 0.66. 


c). Pressure boundary condition 

The static pressure p along the inlet face can be prescribed. The rest of the 
flow variables are extrapolated from the upstream flow field. 

4). Nacelle exhaust plane/propeller disk boundary conditions 

The nacelle exhaust plane as well as the downstream side of the propeller disk 
can have similar boundary conditions. In either case, a subsonic inflow boundary 
is assumed. 

4(a) The first type of boundary condition that can be applied to both nacelle 
exhaust and propeller disk prescribes the total pressure, total temperature and the 
swirl on boundary. These quantities are then used to evaluate the Mach number, 
static pressure and density according to the following formulae 


M = [2q 2 /(y-l)(2Y(T„) 2 /(Y-l)-q 2 )] 1/2 , 

P = p 0 /[l + (Y-l)M 2 /2]^'> , (C6) 
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p = YP/[(Y-l)(C p T 0 - q 2 /2)] 


C p is the specific heat at constant pressure. The subscript o refers to "total" 
quantities. The flow speed q is obtained by extrapolation from the latest iterative 
value from downstream. The swirl distributions are then used to give the u, v, and 
w velocity components from q. 

4(b) Another type of boundary condition frequently used for propeller 
simulation prescribes the distribution of the thrust per unit area (F x ), the normal 
force per unit area (F y ), the side force per unit area (F z ), and the work done by the 
propeller (Q) to calculate the flow variables downstream of the propeller disk. 
These quantities along with the velocity extrapolated from the flow field are used 
in the momentum equations to solve for p, u, v, w and p, downstream of the disk. 


pUu + S x p = P 1 U 1 U 1 + S x pi + SF X , 
pUv + S y p = PiUjVi + SyPi + SFy , 

pUw + S z p = P 1 U 1 W] + S z pi + SF Z , (C7) 

u 2 + v 2 + w 2 = q 2 , 

P = YP/I(Y-1XC p T 0 , + SQ/pU - q 2 /2)] 


where U denotes the contravariant velocity component in the ^-direction with 
subscript 1 denoting upstream conditions. S x , S y , and S z are the projected areas in 
the Cartesian directions of a cell face on the propeller disk. . S is the cell face area. 

In either type of propeller disk boundary condition method, the boundary 
condition on the upstream side of the propeller disk is obtained by setting the mass 
flux to be equal to the mass flux on the downstream side. Here, the subsonic 
inflow assumption implies that both the upstream Mach number Mi and the 
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downstream Mach number M 2 are less than unity. In practice (Ref. 30) transient 
supersonic flow could sometimes appear at the upstream side of the disk during the 
early cycles of the iteration process. When Mj is greater than one, no boundary 
condition can be assigned along the upstream side of the disk. The flow variables 
needed for the calculations are obtained from upstream extrapolation. The 
resulting normal mass flux is then used as a boundary condition for the 
downstream side of the disk. For a subsonic downstream flow (M 2 < 1) this 
implies that one boundary condition described in either type of method must be 
deleted in order to have a well posed boundary-value problem. In the case of 
Method 4(a), the total temperature T 0 is deleted, whereas in the case of Method 
4(b), the work term Q is deleted. In this way, the primary influence of the 
propeller (the thrust) can still be properly simulated. For test cases that are 
representative of candidate installations, such transient supersonic flow no longer 
exists during the latter stage of iteration process, and the proper boundary 
conditions discussed in Methods 4(a) and 4(b) are reinstated. The present 
propeller theory is not valid for the case where both Mj and M 2 are greater than 
one. 
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Fig. 2. Surface Grid for an Engine-Pylon-Airframe Configuration. 
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a. Collapsed Edge b Fictitious Comers 

Fig. 3. Singular or Irregular Grids for a General Multiblock Solver. 



Fig. 4. Surface Grid for a Wing-Mounted Propfan 
in a Wind Tunnel. 
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5a. Configuration. 


5b. Nacelle slices. 



5c. Grid slice in the y-z plane. 


5d. Detail in wing/nacelle region. 



sill 



5e. Grid on nacelle symmetry plane. 


5f. Grid surface at propeller tip. 


Fig. 5. Configuration and Grid Slices. 
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Fig. 7. Propeller Power Effects on Surface C p , Moo=0.167, oc=0.0 deg. 
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Fig. 8. Comparison of Experimental and Computed C p Distributions for 
Propeller Power-On, = 0.167, a = 0.0 deg. 
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Fig. 10. GMBE computed Isobars for Flow Through Nacelle, Moo = 0.77, 
a = 0.50 deg. 
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